clear all;
%% Code explanation
% 170305 DZ 
% 190709 DZ updated to calculate 50%Rn and 1%Rn

%% Define the resistance to temperature conversion equation
T_RuOx=@(R)1./exp(-1.001200439904+1.175279377873*log(R-2.2)+0.021010648114*(log(R-2.2)).^2-0.024007137972*(log(R-2.2)).^3+0.002054398372*(log(R-2.2)).^4);
%
%% Set the parameter

Boffset=0.008;
L_W=24.2/53.3;

ratio=0.5;
%% Read the data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% T=75;
%  T=[10 25 37 50 62 75 87 100 112 125 137 150 175 200]';

root=('D:\Documents2\Rawdata\Sn\171219 MPI mK\20180330 2Sn15PbTe\'); %

filenames={'20180315_0uW_wroots';'20180314_20uW_wroots';'20180314_40uW_wroots';...
    '20180313_60uW_wroots';'20180313_80uW_wroots';'20180312_100uW_wroots';...
    '20180311_130uW_wroots';'20180310_160uW_wroots';...
    '20180309_200uW_wroots';'20180309_250uW_wroots'; '20180316_50uW_wrotary';'20180308_300uW_wroots';...
    '20180317_100uW_wrotary';'20180317_150uW_wrotary';'20180320_400uW_wroots';...
    '20180318_220uW_wrotary';'20180318_300uW_wrotary';'20180321_550uW_wroots';'20180320_650uW_wroots';'20180319_650uW_wroots';'20180319_550uW_wrotary';...
    '20180320_900uW_wroots';'20180319_650uW_wrotary';...
    '20180320_1100uW_wroots';'20180321_1300uW_wroots'};

NT=length(filenames);


colorset=jet(NT);
figure

Tmc=zeros(NT,1);
Tmc_err=zeros(NT,1);
Hc2=zeros(NT,1);
Rc=zeros(NT,1);

%% obtain normal state resistance
kk=NT;
 file1=[root,filenames{kk},'_outofplane_measdown.dat'];
        fid=fopen(file1);
        data1 = textscan(fid,'%f %f %f %f','commentStyle','#');
        data1=data1';
       fclose(fid);
    
    
          [Bmin,i1]=min(data1{1});
          B1=data1{1}(1:i1);
          r1=data1{2}(1:i1)/L_W/50;        
          [Bmin1,ii1]=min(abs(B1));
          B=B1(1:ii1)-Boffset;
          R1=r1(1:ii1);
          
          indexB=find(B>0.5);              
          p1=polyfit(B(indexB),R1(indexB),1);

     plot(B,ratio*(p1(1)*B+p1(2)),'-r');
                  hold on     
          
%%           
for kk=1:NT
   
    file1=[root,filenames{kk},'_outofplane_measdown.dat'];
        fid=fopen(file1);
        data1 = textscan(fid,'%f %f %f %f','commentStyle','#');
        data1=data1';
        fclose(fid);
    
    
          [Bmin,i1]=min(data1{1});
          B1=data1{1}(1:i1);
          r1=data1{2}(1:i1)/L_W/50;
          T=T_RuOx(data1{3}(1:i1)*100);
          Tmc(kk)=mean(T);
          Tmc_err(kk)=3*std(T);
          
          [Bmin1,ii1]=min(abs(B1));
          B=B1(1:ii1)-Boffset;
          R1=r1(1:ii1);
          T=T(1:ii1);

          if kk==17
                 plot(B,R1,'Color',colorset(kk,:)','LineWidth',1);
                  hold on
          end
          if (kk~=11)||(kk~=13)||(kk~=19)
              if (mod(kk,2)==0)&&(kk~=18)
                  plot(B,R1,'Color',colorset(kk,:)','LineWidth',1);
                  hold on
              
              elseif (kk==1)||(kk==25)
                       plot(B,R1,'Color',colorset(kk,:)','LineWidth',1);
                       hold on
              end
          end
%           subplot(1,3,2),plot(data1{1}(index),data1{3}(index),'Color',colorset(kk,:)','LineWidth',ii);
%           hold on
%           subplot(1,3,3),plot(data1{1}(index),data1{4}(index),'Color',colorset(kk,:)','LineWidth',ii);
%           hold on
          
%          indexB=find(B>0.5);
%           p1=polyfit(B(indexB),R1(indexB),1);
         index1=find((R1-ratio*p1(1)*B-ratio*p1(2)<=0));
         if isempty(index1)==0
              Hc2(kk)=B(index1(1));
             Rc(kk)=R1(index1(1));
         else
             Hc2(kk)=NaN;
             Rc(kk)=NaN;
         end
      title={'B(T)','R(kOhm)','T(K)','averaged T(K)'};
        AA=[B R1 T Tmc(kk)*ones(length(B),1)];
       
         xlswrite('FigS6A.xlsx',AA,kk);
       xlswrite('FigS6A.xlsx',title,kk);
%     
end


pp1=polyfit(Tmc(20:25),Hc2(20:25),1);
Tmc1=abs(pp1(2)/pp1(1));

xi=sqrt(6.626E-34/4/pi/1.6E-19/(abs(pp1(2))))*1E9;

Tp=(0:0.002:Tmc1)';

figure
plot([Tmc' Tmc1]',[Hc2' 0]','or',Tp,pp1(1)*Tp+pp1(2),'-k');
axis([0 0.4 0 0.08]);
hold on






%% End




